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ABSTRACT 

The stellar dynamics of u Centauri are inferred from the radial velocities 
of 469 stars measured with CORAVEL ( Mayor et al. 19971) . Rather than 
fit the data to a family of models, we generate estimates of all dynamical 
functions nonparametrically, by direct operation on the data. The cluster is 
assumed to be oblate and edge-on but mass is not assumed to follow light. 
The mean motions are consistent with axisymmetry but the rotation is not 
cylindrical. The peak rotational velocity is 7.9 km s _1 at ~ 11 pc from the 
center. The apparent rotation of u Centauri is attributable in part to its proper 
motion. We reconstruct the stellar velocity ellipsoid as a function of position, 
assuming isotropy in the meridional plane. We find no significant evidence for 
a difference between the velocity dispersions parallel and perpendicular to the 
meridional plane. The mass distribution inferred from the kinematics is slightly 
more extended than, though not strongly inconsistent with, the luminosity 
distribution. We also derive the two-integral distribution function f(E,L z ) 
implied by the velocity data. 
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1. Introduction 

Since the mid-1960's, the standard method for analyzing globular cluster data has been 
a "model-building" approach ( |King 195^| ; |Gunn fc Griffin 1979| ). One begins by postulating 



a functional form for the phase-space density / and the gravitational potential $; often the 
two are linked via Poisson's equation, i.e. the stars described by / are assumed to contain 
all of the mass that contributes to $. This / is then projected into observable space and 
its predictions compared with the data. If the discrepancies are significant, the model is 
rejected and another one is tried. If no combination of functions {/, $} from the adopted 
family can be found that reproduce the data, one typically adds extra degrees of freedom 
until the fit is satisfactory. For instance, / may be allowed to depend on a larger number 
of orbital integrals ( [Lupton, Gunn fc Griffin 1987| ) or the range of possible potentials may 
be increased by postulating additional populations of unseen stars ( pa Gosta fc Freeman 

mm. 



This approach has enjoyed considerable popularity, in part because it is computationally 
straightforward but also because (as |King 198T[ has emphasized) globular cluster data are 



generally well fit by the standard models. But one never knows which of the assumptions 
underlying the models are adhered to by the real system and which are not. For instance, a 
deviation between the surface density profile of a globular cluster and the profile predicted 
by an isotropic model is sometimes taken as evidence that the real cluster is anisotropic. But 
it is equally possible that the adopted form for f(E) is simply in error, since by adjusting 
the dependence of / on E one can reproduce any density profile without anisotropy. Even 
including the additional constraint of a measured velocity dispersion profile does not greatly 
improve matters since it is always possible to trade off the mass distribution with the 
velocity anisotropy in such a way as to leave the observed dispersions unchanged (|Dejonghe 



fc Merritt 1992 ). Conclusions drawn from the model-building studies are hence very difficult 



to interpret; they are valid only to the extent that the assumed functional forms for / and 
$ are correct. 

These arguments suggest that it might be profitable to interpret kinematical data 
from globular clusters in an entirely different manner, placing much stronger demands on 
the data and making fewer ad hoc assumptions about / and $. Ideally, the unknown 
functions should be generated nonparametrically from the data. Such an approach has 
rarely been tried in the past because of the inherent instability of the deprojection process. 
But techniques for dealing with ill-conditioned inverse problems are now well developed 
(e.g. |Wahba 1990| ; [Scott 199"!% |Green fc Silverman 1994p , and the only thing standing in the 



way of a fully nonparametric recovery of / and <£> is the limited size of most kinematical 
data sets. 
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As we show here, the 469 stellar velocities measured by CORAVEL in u 
Centauri ( [Mayor et al. 1997] , Paper I) permit the dynamical inverse problem for 
this system to be solved with relatively few ad hoc assumptions, u Centauri is a 
particularly good candidate for such a study since the two-body relaxation time tft is of 
order 10 9 years even in its core, and tR exceeds a Hubble time at the half-light radius (e.g. 
Meylan et al. 19951 ). Thus there is little justification for assuming that the dynamical state 
of uj Centauri is well described by a model derived from the theory of stellar encounters. 
Instead, we expect the structure of this cluster to reflect in large measure the details of its 
formation process, about which very little is known. The inverse-problem approach is also 
well suited to uj Centauri since the CORAVEL data extend, more or less uniformly, over 
the entire image of the cluster (see Figure 1). 

Our analysis will proceed by discrete steps, each of which takes us deeper into the 
dynamics though at the price of a progressively larger set of assumptions. As a first step 
we recover the rotational velocity field, assuming only that uj Centauri is axisymmetric 
and edge-on. (The assumption that we lie in the equatorial plane of uj Centauri is made 
throughout this paper, partly for reasons of mathematical convenience, but also because 
uj Centauri is among the flattest of globular clusters and hence it is unlikely that we 
are viewing it from a direction far out of the equatorial plane.) We next recover the two 
components cx ro and of the stellar velocity ellipsoid, after imposing dynamical equilibrium 
and assuming isotropy in the meridional plane, i.e. a m = a z . The latter assumption is 
the least restrictive one we can make which still allows us to uniquely recover the internal 
velocity dispersions from the observed ones. The two-dimensional gravitational potential 
and the mass density then follow without additional assumptions; we emphasize that these 
are bona-fide dynamical estimates and are independent of any assumptions about the 
mass-to-light ratio or the stellar mass function. Finally, we recover / by assuming that 
the phase space density depends only on the orbital energy E and angular momentum L z 
- consistent with our previous assumption that a w = <r z , but slightly more restrictive in 
the sense that even three-integral /'s can be isotropic in the meridional plane. At no point 
do we write down ad hoc forms for any of the unknown functions; the only restriction we 
impose is smoothness, enforced via penalty functions. Details of the technique may be 
found below and in Merritt (1996). 

One advantage of analyzing the data in this way is that the conclusions drawn at every 
step are much easier to interpret than if / and $ were written down at the outset. For 
instance, the dependence of the mean streaming velocity on position in the meridional 
plane follows immediately from the velocity data without any assumptions other than 
symmetry, and the uncertainties in our recovery of this function are due almost entirely 
to finite-sample fluctuations (which we estimate via the bootstrap) - there are no biases 
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like those that would result from the adoption of a rigid functional form for or /. In 
the same way, our conclusions about the degree of velocity anisotropy in u Centauri are 
based on direct reconstruction of the velocity ellipsoid, and not merely on deviations of the 
density profile from the predictions of an ad hoc isotropic model. 

By proceeding in this step-by-step way, we also preserve the important distinction 
between the "knowability" of different dynamical quantities. For instance, the phase-space 
density / is related to the data by what is effectively a high-order differentiation, which 
means that recovery of / from kinematical data will always be extremely ill-conditioned 
( [Miller 1963| ). It is therefore proper to save the determination of / until the last, after 
quantities that are more directly related to the data - like the rotational velocity field 
- have been derived. This distinction between well-determined and poorly-determined 
quantities is lost in the model-building approach, which begins by specifying / and then 
treats all quantities that are derived from it on an equal footing. 

We find that to Centauri can be well described as an "isotropic oblate rotator," i.e. a 
system in which the velocity residuals about the mean motions are approximately isotropic. 
Our conclusion is stronger than one derived from the tensor virial theorem alone ( [Meylan 



fc Mayor 1986|), since we are able to make statements about the detailed dependence of 



the velocity ellipsoid on position. We find in addition that the mass distribution in u 
Centauri is consistent with the luminosity distribution of the bright stars, and that the 
mass density normalization is consistent with that predicted by the virial theorem under 
the assumption that mass follows light. Hence there is no compelling reason to postulate 
a population of unseen objects in this cluster. Finally, we show that the velocity data are 
fully consistent with a two-integral stellar distribution function f(E,L z ). 

In §2 we review the CORAVEL data, which are described more fully in Paper I, and 
briefly describe the penalized likelihood formalism that is used in the reconstruction of the 
dynamical quantities. The stellar density profile is derived in §3; the rotational velocity 
field in §4; and the velocity dispersions in §5. Dynamical estimates of the gravitational 
potential and the mass distribution are presented in §6, and the stellar distribution function 
is derived in §7. §8 sums up and discusses prospects for future work. 



2. Data and Method 



The data on which this study is based are presented in Paper I. The photoelectric 
spectrometer CORAVEL ([Mayor 1985|) was used to obtain 1701 radial velocity measurements 
of 483 bright stars, mostly giants and subgiants, of which 469 were later determined to 
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be members. The typical measurement uncertainty is < 1.0 km s _1 . The positions of the 
member stars are shown in Figure 1. Previous analyses of various subsets of these data 
have been presented by Meylan & Mayor (1986), Meylan (1987) and Meylan et al. (1995). 

We adopt Cartesian coordinates (X, Z) on the plane of the sky; Z is coincident with 
the isophotal minor axis of uj Centauri, as defined below, and X increases westward. The 
coordinate parallel to the line of sight is then Y. We take minutes of arc as our units for X 
and Z. Adopting a heliocentric distance to to Centauri of 5.2 kpc ( |Meylan 1987| ), we find 



1' = 1.51pc. The canonical "core radius" of uj Centauri is ~ 3' ~ 4 pc ( Peterson fc King 



|1975| ) (although we show below that the brightness profile of uj Centauri is consistent with 
a luminosity density that increases monotonically toward the center). The tidal radius r t is 
roughly 50' ~ 75 pc ( [Peterson fc King 1975| ). For comparison, the velocity data in Paper I 
extend to about l/2r t , or ~ 25' (Figure 1). 

We assume throughout most of this study that uj Centauri is axisymmetric. Cylindrical 
coordinates in the meridional plane are (zcr, z); the z axis is parallel to the Z axis, i.e. 
coincident with the isophotal minor axis. 

Our treatment of the data will follow the nonparametric approach developed in an 
earlier series of papers ( [Merritt 1993b| , 1996; Merritt fc Tremblay 1994] ). One seeks a 



smooth function - call it g - such that the projection of g into observable space is consistent 
with the data. For instance, g might be the rotational velocity field in the meridional plane, 
to be determined from a set of measured positions and radial velocities. The standard 
technique for solving such inverse problems is to vary g so as to minimize a functional like 



- log C p = 



data 



9 P -Ag' 1 



+ aP(g), (1) 



the "penalized log likelihood" ( Thompson fc Tapia 1990[) . Here A represents the projection 



operator that brings g into observable space; g p is the observed quantity; e is a measurement 
error; and P is a function that assigns a large penalty to noisy solutions. The penalty 
function is needed since direct deprojection of the data is generally ill-conditioned, leading 
to physically useless solutions ( [King 198T| ) . In the implementations presented here, P will 
depend on g via its mean square second derivatives, a fairly standard choice (e.g. |Wahba| 



|fc Wendelberger 1980| ). The degree of smoothness of the solution is then controlled by 
adjusting a. One generally chooses a to be as small as possible consistent with smoothness 
so as to avoid biasing the solution. 

Implicit in this approach is the assumption that the observable quantities g p contain 
sufficient information to uniquely constrain the unknown functions g. The requisite proofs 
for the axisymmetric inverse problem are presented in Merritt (1996). 
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3. Luminosity Density Profile 

The first inverse problem to be solved is the derivation of the spatial luminosity 
density in u Centauri, u(zj,z), as a function of position in the meridional plane, given 
measurements of the surface brightness £(JT, Z) on the plane of the sky. Ideally we would 
adopt for £ the density of stars that comprise our kinematical sample, computed via a 
penalized-likelihood scheme flMerritt fc Tremblay 19941) . However the velocity sample was 
not specifically chosen to be magnitude-limited (see Paper I), nor is it as large as one would 
like for the nonparametric reconstruction of v. Instead we follow the practice of Meylan 
(1987) and use surface brightness determinations from a number of other sources. These 
include centered aperture photometry ( Gascoigne fc Burr 1956 ; Da Costa 1979|) and drift 



scan measurements ( pa Costa 1979| ) in the central regions, and star counts ( [King et al. 



1968[ ) in the outer regions. The normalized surface brightness measurements and their 
estimated errors are presented in Table 1 of Meylan (1987) at 44 distinct radii i?,. The 
studies from which these measurements were taken generally ignored the noncircular shape 
of L) Centauri, and so the radii tabulated by Meylan (1987) and used here should be viewed 
as approximate averages over position angle. 

uj Centauri is however significantly flattened. Geyer, Hopp & Nelles (1983) measured 
the isophotal ellipticity e = 1 — b/a as a function of radius between 1.4' and 28' via 
photographic photometry. They found a mean ellipticity of 0.121; however e varies 
significantly with radius, from e ~ within 2' to e ~ 0.25 at 10', becoming rounder again at 
large radii. The position angle of the principal axis was not found to vary significantly with 
radius. Geyer et al. do not quote a value for this angle, but White & Shawl (1987) derive 
an orientation of 6° east-from-north for the isophotal minor axis, and this value appears to 
be consistent with Figure 1 from Geyer et al. We will adopt the White & Shawl orientation 
in what follows. 

Ideally, one would carry out a full deprojection of the surface brightness data £j, 
measured at some set of points {X i: Zj}, to obtain an estimate of the space density u(w, z). 
This inverse problem has a formally unique solution for any axisymmetric system that is 
viewed edge-on, even if its cross-section is not elliptical ( [Rybicki 1986j ; |Gerhard fc Binney 



|1996| ). However no two-dimensional surface brightness data for u Centauri have been 
published. Instead, we will assume that the density of stars in uo Centauri is stratified on 
similar, concentric oblate spheroids of axis ratio b/a = 0.879, the average value found by 
Geyer et al. (1983) — consistent with our assumption that the equatorial plane is parallel 
to the line of sight. 

The binned surface brightness measurements from Table 1 of Meylan (1987) are plotted 
in Figure 2a. The solid line in that figure is an estimate of the surface brightness profile, 
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defined as the function E(i?) that minimizes the quantity 
- log C p = £ er 2 [log S< - log S(^)] 2 + a | 



^ 2 (logS) 
c/(lo gj R) 2 



dlogi?. (2) 



Here is the estimated uncertainty in the surface brightness measurement logSj at i?,. 
The value of a used to derive the profile of Figure 2a was selected via the generalized 
cross-validation technique flWahba 1990| , ch. 4); this value minimizes the estimated, 
integrated square error in S(i2). 

The estimate of v{r) may be defined as the Abel inversion of the estimate S(i£): 

* , N 1 f°° dR . . 

Here r is also an azimuthally-averaged mean radius. But deprojection has the property of 
amplifying the noise in the data, and we expect that the optimal value of a to be applied 
to the surface brightness data when deriving v will be larger than the optimal value for 



estimating E itself (e.g. |Scott 1993 , p. 132). Unfortunately, there is no generally-accepted 
criterion for selecting a in cases like this ( [Wahba 1990| , p. 105). To produce the estimate 
of u(r) shown in Figure 2b, we used a value of a roughly three times the optimal value for 
estimation of S(i2). The 95% confidence bands in that figure were derived via the bootstrap 
( |Wahba 1990| , p. 71), with a pointwise correction for bias (as in |Scott 199~2| p. 259). 



The stellar density profile u(r) is evidently very well determined by these data at 
all radii outside of ~ V. (For comparison, Peterson fc King 1975| find ~ 3' for the u 
Centauri core radius.) However at smaller radii the density profile is poorly determined, 
due primarily to the amplification of uncertainties resulting from the deprojection. The 
profile of Figure 2b actually has a power-law cusp, v ~ r^ 1 , inside of 0.5'; however the 
confidence bands are consistent with a wide range of slopes in this region, including even a 
profile that declines toward the center. 

Henceforth we take as our estimate of the stellar number density v{w,z) = with 
£ 2 = (1 - e) [zu 2 + z 2 /(l - e) 2 ] and v the function of Fi Kure 2b. Almost all of the stars with 
measured velocities lie outside of the region where the uncertainties in v are significant; 
furthermore, the uncertainties in most of the quantities derived below will be affected much 
more by limitations in the kinematical sample than by errors in v. Hence we will ignore 
uncertainties in v in what follows. 
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4. Rotational Velocity Field 

Next we wish to find the dependence of the mean azimuthal velocity on position in 
the meridional plane. As a first step, we investigate the variation of the mean line-of-sight 
velocity Vy(X, Z) over the image of uj Centauri. An estimate of Vy can be defined as the 
function Vy that minimizes 

2 



log£ p = £e7 2 V t -V Y {X h Z, 



+ a 



'dWy\ 2 | 2 (d 2 Vy\ 2 | (dW, 



dX 2 / \dXdZ \ dZ 2 



dXdZ. (4) 



Here ej is the estimated uncertainty in the stellar velocity V, measured at point (Xj,Zj). 
The function Vy that minimizes (f|) is a so-called "thin-plate smoothing spline" ( |Wahba| 



1990| , p. 30). Figure 3a shows the result; the smoothing parameter a was again chosen via 
generalized cross-validation. The rotational velocity field is approximately symmetric but 
shows what appears to be a systematic twist at large radii. 

Here we note a correction that must be made to the measured velocities before 
proceeding further, uj Centauri has a substantial proper motion. Feast, Thackeray & 
Wesselink (1961) note that the translation of a solid object produces an apparent rotation, 
since the projection of the space velocity along the line of sight is different at different 
points in the image. This "perspective rotation" results in a radial velocity increment 

V pr = ^V x + |Vz (5) 

where Vx and Vz are the components of the cluster space velocity parallel to the X and Z 
axes and D is the distance to the cluster. Curves of constant V pr are straight lines of slope 
—Vx/Vz, and the apparent "rotation axis" is oriented at an angle 9 pr = tan _1 (— Vx/Vz) 
measured clockwise from the Z-axis. As a result of the perspective rotation, the cluster 
appears to rotate as a solid body, with the rotational velocity increasing linearly with 
radius. 

The proper motion of uj Centauri has been determined by Murray, Jones & Candy 
(1965). Corrected for differential Galactic rotation (Cudworth 1994), the components of 
the proper motion parallel to our X and Z axes are 

fix = 0.26 ± 0.06 arcsec/century, fi Z = —0.74 ± 0.05 arcsec/century. (6) 

Using our adopted distance of 5.2 kpc and assuming an error in this number of ± 10%, we 
find 

V x = 64.1 ± 16 km s~\ V z = -183.0 ± 22 km s~\ (7) 
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The resulting change in the radial velocity at point (X, Z) is 

(0.0187 ± 0.005 X' - 0.053 ± 0.006 Z') km s" 1 (8) 

with X' and Z 1 measured in minutes of arc from the cluster center. The amplitude of the 
effect is small but not negligible, exceeding 1 km s _1 near the edge of our kinematical 
sample. The spurious "rotation axis" is oriented at 19 degrees west-from-north, or about 25 
degrees in a clockwise direction from the isophotal minor axis. At points near the X axis, 
the sign of the induced radial velocity is opposite to that of the true, line-of-sight rotational 
velocity. Thus the perspective rotation would be expected to induce a slight twist at large 
radii in the contours of constant mean radial velocity — and in just the direction observed 
in Figure 3a. 

Figure 4a shows the corrected rotational velocity field, computed by minimizing (f|) 
after removing the contribution of the perspective rotation to the measured velocities. The 
contours are now substantially more symmetric at large radii. The rotational velocity along 
the isophotal major axis, shown in Figure 4c, has a peak value of about 7 km s _1 at a 
radius of 7' ~ 11 pc; this profile is reasonably symmetric about the cluster minor axis. 

If uj Centauri were exactly axisymmetric, the kinematic minor axis would coincide 
with the isophotal minor axis. Figure 4a shows that this is approximately true, although 
there still appears to be a separation of a few degrees between the axis of maximal apparent 
rotation and the isophotal major axis (the X-axis in Figure 4). The discrepancy is probably 
not significant given the relatively small number of velocities that determine our estimate 
of the rotational velocity field, and given the likely error in White & Shawl's (1987) 
determination of the minor axis. Furthermore the proper motion adopted here might be in 
error, leading to an incorrect adjustment for the perspective rotation. In any case, we will 
continue to use the White & Shawl estimate of the minor axis orientation in what follows. 

Our next task is to recover v<p(m,z), the mean rotational velocity in the meridional 
plane. This inverse problem has a formally unique solution if uj Centauri is axially 
symmetric and edge-on ( [Mcrrificld 1991 ). However because is related to the data via a 



deprojection, it will be intrinsically less well determined than Vy. 
We seek the function v^w, z) that minimizes 

- log C p = £i 2 {Vi ~ £iV<p) 2 + a 




dw 2 i - i ;) _ ;) . i i :>. > i dwdz, ( 9 ) 
where £ is the line-of-sight projection operator 

u(w, Z^v^tu, Zi) —j= (10) 



-10- 



Equation (|T0| ) assumes that the cluster is axially symmetric, and that the observer lies in 
the equatorial plane. We can decrease the noise in the estimate by assuming in addition 
that the rotational velocity field is symmetric about the equatorial plane, an assumption 
that makes dynamical sense and that is consistent with the appearance of the contours in 
Figure 4. One way to enforce this constraint is to reflect all of the data into one quadrant 
of the (X, Z) plane, with appropriate changes in the sign of V^. This procedure was found 
to give reasonable results, except that the estimate so produced had a discontinuous 
gradient along the w and z axes. Physically more reasonable results were obtained by 
changing only the signs of the X i; and requiring (via a set of linear constraints) that the 
solution be symmetric about z = and zero along the z axis. The penalty function was 
then able to enforce nearly-zero ^-derivatives on the equatorial plane, as desired. 

The corresponding, constrained optimization problem was solved via quadratic 
programming flMerritt 1996 ). Various values were tried for the smoothing parameter; 



Figure 5 shows the result for one choice of a judged nearly optimal. The contours of 
constant are remarkably similar in shape to those of the parametric model postulated 
by Meylan & Mayor (1986), at least in the region near the center where the solution is 
strongly constrained by the data. As noted by them, the rotational velocity field is clearly 
not cylindrical; instead, has a peak value of ~ 8 km s" 1 at about 7' from the center in 
the equatorial plane, and falls off both with increasing w and z. In the region inside the 
peak, the rotation is approximately solid-body. 

One would like to estimate the uncertainties associated with the solution presented in 
Figure 5. The error in can be broken into two parts: the "variance," i.e. the mean square 
fluctuation resulting from the finite sample size and measurement errors; and the "bias," 
the systematic deviation of the solution from the true (e.g. jMiillc-r 1988| , p. 29). The 



variance is possible to estimate via the bootstrap, by generating a large number of pseudo 
data sets from the smooth solution and observing how greatly the estimates of fluctuate 
from sample to sample. The bias is a tougher nut to crack. While the techniques used 
here are (unlike parametric methods) unbiased in the limit of large sample sizes, for finite 
samples the smoothing will introduce a systematic deviation of the solution from the true 
Vfj,. The situation is complicated still more by the fact that the amplitude of the variance 
depends (inversely) on the amount of smoothing, and hence on the bias. Thus a too-large 
choice for a will artificially reduce the variance that would otherwise be expected from a 
sample of a given size. (Exactly the same is true for a histogram - the y/N uncertainty 
associated with any bin depends on the bin width, i.e. on the degree of smoothing.) 

We nevertheless constructed formal confidence bands on our estimate of according 
to the following scheme ([Wahba 1990| , p. 71). We emphasize that intervals constructed in 
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this way are confidence intervals for the estimate and not for the true function v^, and 
do not reflect the error due to the bias. We expect the latter to be greatest at large radii, 
where the data constrain the solution the least. 

1. Compute a smooth estimate of v^; in our case, this estimate is the function plotted in 
Figure 5. 

2. Project this estimate into observable space and record the predicted, mean velocity V \ 
at the positions (X i} Z/) of the observed stars. 

3. Generate a large number N of pseudo data sets Vi = Vi + ii, where q is a random 
number from the normal distribution A/"(0, of), and of is an estimate of the variance of the 
measured velocities about their mean value at (Xj, Z^). This variance was set equal to the 
square of the estimated velocity dispersion displayed in Figure 4b. (Measurement errors in 
the Vi are of order 1 km s _1 and are negligible compared to the intrinsic velocity dispersion.) 

4. For each of these pseudo data sets, compute a new estimate of v<f> and record the values 
at each point on the solution grid. 

5. The 5 confidence intervals at any point on the grid are given by the 5/2 Nth. and 
(1 — 5/2) iVth values from the distribution of estimates v<p at that point. 

6. Since the pseudo data are generated from the smoothed estimate v^, which is biased, the 
estimates generated from the pseudo data will be "doubly" biased. One can compute the 
bias between the original estimate and the bootstrap estimates and then subtract it from 
the computed confidence bands (e.g. [Scott 1952] , p. 259). This (slightly ad hoc) procedure 
yields confidence bands that are correctly situated with respect to the original estimate, 
rather than with respect to the average of the bootstrap estimates. 

Figure 6 gives 95% confidence intervals so constructed on v^, along two lines in the 
meridional plane: as a function of w at z = 0, and as a function of z at the value of w 
corresponding to peak rotation in the equatorial plane. The confidence intervals on are 
wide, partly a consequence of the fact that relatively little of the motion in oj Centauri is 
ordered. The peak rotational velocity is 7.9^2i5 km s" 1 (95%). For comparison, Meylan & 
Mayor (1986), using an ad hoc parametric form for v^, derived a peak value of ~ 8 km s _1 
at ~ 7', in excellent agreement with our result. Outside ~ 15', the confidence bands shown 
in Figure 6 suggest that the form of the rotational velocity field is only weakly constrained 
by these data. 

We would like to compare our inferred rotational velocity field in u Centauri to the 
predictions of a theoretical model. However it is not clear that any relevant models exist. 
Because of the long relaxation time in u Centauri, the rotation probably still reflects to a 
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large extent the state of the cluster shortly after its formation. Our estimate of v^(w, z) 
might therefore be most useful as a constraint on cluster formation models. 

A simpler question is whether uj Centauri can be described as an "isotropic oblate 
rotator" ; that is, whether the distribution of velocity residuals about the mean motion are 
approximately the same in all directions. We will answer this question in the affirmative 
below, by showing that there is no evidence for a significant difference between the velocity 
dispersions parallel and perpendicular to the meridional plane. Hence u Centauri can 
reasonably be described as "rotationally flattened." Meylan & Mayor (1986) reached a 
similar conclusion after noting a strong correlation of the local rotational velocity with the 
isophotal flattening. 



5. Velocity Dispersions 



We now wish to dig deeper into the dynamics and reconstruct the dependence of the 
velocity dispersions on position in the meridional plane. We again start by investigating the 
variation in the plane of the sky. We seek the function Vy(X, Z) that minimizes 



-log£ p = 52e?\v?-V?(Xi,Z, 



+ a 



'd 2 Vf 
dX 2 



2 

\ dXdZ / 



'd 2 Vf 
dZ 2 



dXdZ; 



then <7y — Vy — V Y . The result is shown in Figures 3b and 4b, before and after removal of 
the perspective rotation. The contours of constant ay are reasonably symmetric, at least 
near the center; however there is no clear indication of an elongation of the contours in the 
X-direction, parallel to the isophotal major axis. The velocity dispersion along the major 
axis (Figure 4d) falls from a maximum of ~ 17 km s _1 at the center to ~ 8 km s _1 at 15'. 

Meylan et al. (1995, Figure 1) present a circularly-symmetrized velocity dispersion 
profile derived from the same data used here. They note that the velocity dispersion rises 
sharply very near the center of u Centauri, from ~ 17 km s -1 at 2' to ~ 22 km s" 1 within 
1'. This central value is derived from only 16 stars and has a formal (la) uncertainty of 
±3.9 km s -1 ; thus the increase near the center may not be significant. We can mimic the 
effect of a fine radial grid like the one used in Meylan et al. (1995) by reducing the value 
of the smoothing parameter a. We find that our estimate of the central dispersion rises as 
a is reduced, and we can easily reproduce the Meylan et al. value for ay(0) if a is chosen 
to be sufficiently small. However the velocity dispersion profile so produced is very noisy 
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and almost certainly undersmoothed; the high central dispersion has the appearance of a 
finite-sample fluctuation. We conclude that there is no secure case to be made from the 
current data for a sharply-rising velocity dispersion within the inner minute of arc. 

The line-of-sight velocity dispersions contain contributions from both o~ w and 
a z does not contribute due to our assumption that the observer sits in the equatorial 
plane. If we assume in addition that <r ro = o z = o~(w,z), i.e. that the velocity ellipsoid 
is circular in the meridional plane, it becomes possible to infer both a(m, z) and a<p(-w, z) 
independently from Vy{X, Z) ([Merritt 1996|) . The assumption of isotropy in the meridional 
plane is physically restrictive and one would like to avoid making it. However, allowing 
radial anisotropy adds two unknown functions to be determined from the data: a z (w,z) 
and v zo v z {w, z), which together with a m and define the 3-D shape and orientation of 
the velocity ellipsoid. The dynamical inverse problem under this less-restrictive set of 
assumptions has not been formulated but one would certainly need much more data than 
are currently available in uj Centauri in order to find a unique solution. This issue is 
discussed at greater length in §8. 

Following the scheme in Merritt (1996), we search for the functions a 2 and v\ that 
minimize 

-log£ p = E^^-Ap,^} 

'd 2 a 



+ a 



+ a 




dwdz 



dwdz 



(12) 



subject to the constraints 



dv da 2 dv da 2 v d 



dzu dz dz dzu zu dz 







(13) 



which are satisfied for any axisymmetric system with cr ro = a z = a. The projection operator 
Li becomes 



£i p, v 2 } = 2S(X,, Z^- 1 (°° v{w, Z t 



l---^\a 2 (w,Z i 



TJJ 



dw 



lw 2 - Xf 
(14) 

and <j| = v 2 — v±. This is again a quadratic programming problem. The results are shown in 
Figure 7. The vertical gradients of a and as were forced to be zero along the w— axis. There 
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are obvious differences between the contour shapes for a and a^, though both functions are 
reasonably symmetric. 

We would again like to know which of the details of Figure 7 are securely implied by 
the data and which are due to finite-sample fluctuations. Perhaps the most interesting 
question here is the value of the anisotropy parameter 7 = 1 — cr^/a 2 , the degree to which 
the azimuthal and meridional dispersions differ. Figures 8 and 9 give estimates of a, and 
7 on the major and minor axes, along with 95% bootstrap confidence bands. Although 
is less than a at most radii, the confidence bands are wide and there is no clear indication 
of anisotropy anywhere in the meridional plane. Near the center, could be as small as 
60% of a or as large as 1.4<r (95%). The formal 95% confidence interval on o"(0,0) is 17^6 
km s -1 . 

The simplest interpretation of these results is that uj Centauri is an "isotropic oblate 
rotator," i.e. that the velocity residuals about the mean motion are isotropically distributed 
everywhere. However we emphasize the uncertainties associated with this interpretation. 
We have ruled out radial anisotropies by fiat, and even under this restrictive assumption, 
the confidence bands on a and are quite wide. Clearly, a somewhat larger kinematical 
sample than ours would be needed to place reasonably strong constraints on the degree of 
anisotropy in uj Centauri. 

Meylan (1987) noted that anisotropic models from the Michie-King family could be 
made to fit the surface brightness data for uj Centauri only if the velocity ellipsoid was 
allowed to become significantly elongated outside of ~ 10'. Meylan's result highlights the 
limited range of models from the Michie-King (or any similar) family. In the absence of a 
priori information about the likely functional form of f(E, L), the surface brightness profile 
of a stellar system does not contain any useful information about the degree of velocity 
anisotropy. Particularly in clusters like uj Centauri, where the collisional relaxation time is 
long, conclusions derived from comparison with collisional models should probably not be 
given much weight. 



6. Mass Distribution 

Having obtained estimates of the first and second velocity moments of the stellar 
distribution function, we can now construct dynamical estimates of the gravitational 
potential &(m,z) and the mass density p{w,z). The former follows most directly from the 
Jeans equation relating the vertical gradients in the stellar pressure to the vertical force: 
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or 

r°° - d log V 



&(zu,z)= I a 2 (zu,z) — - — dz — a 2 (w, z). (16) 

J z UZ 

We have again made use of our assumption that a ro = o z = a. The mass density follows 
from Poisson's equation, 



p{w,z) 



AttG 



1 d ( d$\ , <9 2 <l 



w dw \ dw I dz 2 



(17) 



Estimates of $ and p so derived are bona-fide dynamical estimates; they are independent 
of any assumptions about the relative distributions of mass and light, the stellar mass 
function, etc. However these estimates are dependent on our neglect of anisotropy in the 
meridional plane, a point that we return to below. 

Figure 10 gives the bias-corrected estimates of $ and p. Also shown in Figure 10 are 
the potential and mass distribution that u Centauri would have if p were proportional to 
z/, i.e. if mass followed light. The total mass of u Centauri under the latter assumption 
follows immediately from the virial theorem (Appendix; |Meylan fc Mayor 1986| ): it is 



1.45 x IQ a (Vy)M & « 2.86 x 10 6 M Q , with (Vy) the mean square line-of-sight velocity 
expressed in km 2 /s 2 . The dynamically-inferred potential is slightly more elongated at large 
radii, r > 5', than the "potential" generated by the light. The differences between p(w, z) 
and z) are somewhat greater, due in part to the fact that p is a second derivative of 
$. The dependence of both functions on radius along the major axis and their confidence 
bands are shown in Figure 11. The dynamically inferred mass distribution is less centrally 
concentrated than the light distribution. However this difference is hardly significant given 
the width of the confidence bands on p. Furthermore, we would expect the smoothing to 
reduce the degree of central concentration in the inferred mass, and our "bias correction" 
probably does not account fully for this systematic error. Thus we are inclined to interpret 
Figure 11 conservatively, i.e. to conclude that there is no evidence for a significant difference 
between the mass and light distributions in uj Centauri. 

In any case, our formal, 95% confidence interval on the central mass density of to 
Centauri is 2110t5foM Q pc~ 3 . This is consistent with an estimate based on the "core-fitting" 
formula, p(0) « 9a Y {0)/A7rGR 2 c (|King 1966; ). Adopting R c = 3.0' = 4.6 pc and a Y = 17 km 



s -1 , we find p(0) = 2270M Q pc~ 3 . This agreement is to be expected since the assumptions 
underlying the core-fitting formula - velocity anisotropy and a constant M/L - are 
consistent with our analysis of the u Centauri data. 

Although we have failed to falsify the mass-follows-light model for u Centauri, we 
stress that very different mass distributions are certain to be consistent with our limited 
data. By far the greatest uncertainty in estimates of $ and p for hot stellar systems results 
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from lack of knowledge about the radial elongation of the velocity ellipsoid. In a spherical 
system, allowing a r ja^ to have any physically permissible value leads to order-of-magnitude 
uncertainties in the total mass or the central density as derived from the velocity dispersion 
profile ( Dejonghe fc Merritt 199^ ). The same is likely to be true in axisymmetric systems, 
with the role of the radial anisotropy played by the elongation of the velocity ellipsoid 
in the meridional plane. If the velocity ellipsoid became radially elongated at large radii, 
for instance, the inferred mass density profile would be less centrally concentrated than 
what we found assuming a w = o z . Possible routes for overcoming this indeterminacy are 
discussed in 58. 



7. Distribution Function 

Our final step is to estimate the formally unique, two-integral distribution function 
f(E, L z ) that generates the inferred stellar number density u(zu, z) and rotational velocity 
field v^{w,z) in the inferred potential &(vj,z). We are justified in expressing / as a 
function of only two integrals since any f(E, L z ) implies a velocity ellipsoid that is isotropic 
in the meridional plane, consistent with our assumption that a w = o z . However there also 
exist three-integral f's with this property and so the distribution function that we derive 
here is not strictly unique even given our restrictive assumption about the kinematics. Our 
goal is simply to show that a non-negative, two-integral / exists that is consistent with our 
estimates of u, and $. The precise functional form of f(E, L z ) is of less interest, for two 
reasons. First, we have already derived the lowest-order moments of /, i.e. v,v^<r and 
Thus / itself tells us relatively little that we do not already know. Second, like the mass 
density p, f is related to the data via a high-order differentiation and hence its derivation 
is extremely ill-conditioned. We expect that many different expressions for f(E, L z ) - 
obtained with different smoothing parameters, or different penalty functions, etc. - will be 
almost equivalent in terms of their ability to reproduce the input functions v and v$ and in 
fact we found this to be the case. 

Figure 12 nevertheless shows estimates of the even and odd parts of /: 

f(E, L z ) = f + (E, L z ) + /_(£, L z ), f ± (E, L z ) = ~ [f(E, L z ) ± f(E, -L z )\ , (18) 



the former derived from u and the latter from ( [Merritt 1996|) . These estimates were 



constrained to be everywhere positive. The positivity constraint did not preclude the 
inferred f's from reproducing the input functions with high accuracy, and we infer from 
this that u Centauri can in fact be well described by a two- integral /. 
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8. Conclusions and Directions for Future Work 

Our principle conclusions are summarized here. 

1. The rotational velocity field in uj Centauri is consistent with axisymmetry, once a 
correction is made for "perspective rotation" resulting from the cluster's proper motion. 
The rotation is strongly non-cylindrical, with a peak rotation speed of 7. 9*2.t km s -1 (95%) 
at a distance of ~ 11 pc from the cluster center in the equatorial plane. The rotation is 
approximately solid-body at small radii; at large radii, the available data do not strongly 
constrain the form of the rotational velocity field. 

2. By assuming that the residual velocities are isotropic in the meridional plane, 
cr ro = a z = a, we derived the dependence of the two independent velocity dispersions a 
and on position in the meridional plane. The central velocity dispersion parallel to the 
meridional plane is <t(0, 0) = 17l|g km s _1 . There is no evidence for significant anisotropy 
anywhere in uj Centauri. Thus, this cluster can reasonably be described as an "isotropic 
oblate rotator." 

3. The gravitational potential and mass distribution in u Centauri are consistent with 
the predictions of a model in which the mass is distributed in the same way as the bright 
stars. The central mass density is 211Oi5ioM pc~ 3 . However this result may be strongly 
dependent on our assumption that the velocity ellipsoid is isotropic in the meridional plane. 

4. We derive a two-integral distribution function f(E, L z ) for the stars in u 
Centauri and show that it is fully consistent with our data. 

Our method is based on the smallest set of assumptions that permit formally unique 
estimates of the first and second moments of the stellar velocity distribution in an edge-on, 
axisymmetric cluster, given observed values of these moments on the plane of the sky. 
The deprojection of the line-of-sight velocity dispersions can be carried out only if the 
velocity dispersions in the meridional plane are reducible to a single function of position 
a(zu,z), rather than the three functions {<r ro , a z , v^v z } that characterize a fully general 
axisymmetric system. While there are many possible choices for the relation between these 
three functions, setting v w v z = and er ro = o z as we do is the only choice consistent with 
a general two-integral distribution function f(E,L z ). Nevertheless, some of our results 
- particularly the form of the mass distribution - might be strongly dependent on our 
neglect of velocity anisotropy in the meridional plane. Future work on uj Centauri should 
therefore be directed toward obtaining enough kinematical data to independently constrain 
the functions a z and v m v z . 

One route would be to measure velocity dispersions parallel to the plane of the sky from 
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internal proper motions. Proper motion velocity dispersions have in fact been measured 
in a few globular clusters (e.g. M13, |Cudworth fc Monet 1979| ; M3, |Cudworth 1979| ). In 



a spherical system, such data allow one to infer the velocity anisotropy as a function of 
radius QLeonard fc Merritt 1989]), and this technique has been used to constrain the mass 



distribution in M13 QLeonard, Richer fc Fahlman 1992)) . The corresponding existence proof 



has yet to be carried out for axisymmetric systems, but we note that proper motion velocity 
dispersions would contribute two new functions of two variables {ax(X, Z), az(X, Z)}, 
which might provide just enough information to uniquely constrain the two unknown 
functions cr z (zu, z) and v VJ v z (w, z). 

Alternatively, one could try to obtain much larger samples of radial velocities in u 
Centauri, of order 10 4 , which would permit the derivation of the full line-of-sight velocity 
distribution N(X, Z, Vy) at a discrete number of positions in the image of the cluster. In 
fact the radial velocities of approximately 4200 stars in the core of ui Centauri have already 
been measured using the Rutgers Fabry-Perot interferometer (|Gebhardt et al. 1997|) , and 
line-of-sight velocity distributions for these stars have been presented by Merritt (1997). 
An observational program directed toward measuring a similar number of radial velocities 
at larger radii in uo Centauri is currently nearing completion ( [Reijns et al. 1993| ). Here 
again, we lack a rigorous theoretical understanding of precisely what dynamical information 
are contained in the function N(X, Z, Vy). However numerical experiments in the spherical 
geometry flMerritt 1993a| ; |Gerhard 1993| ) suggest that line-of-sight velocity distributions, 



like proper motions, contain considerable information about the velocity anisotropy. 

The kinematical sample analyzed here is large by the standards of just a few years ago, 
yet still small enough that the confidence intervals on many quantities of interest are still 
quite wide. We trust that this situation will be rectified in the very near future. 



We thank C. Pryor for comments on the manuscript. This work was supported by NSF 
grant AST 90-16515 and NASA grant NAG 5-2803 to DM. 



-19- 



APPENDIX 

Here we derive an expression for the total mass of Omega Centauri assuming that the 
mass density is everywhere proportional to luminosity density. 

The scalar virial theorem is 2K + W = 0, with K the total kinetic energy and W the 
potential energy. Since we have assumed that the starlight, and hence the mass, is stratified 
on similar oblate spheroids, W can be written 



W = -An 2 GS ( - ) — — — sin^e (19) 



2 



a 



( [Roberts 1962| ) where 

S = / dm 2 p(m 2 ) / p(m ,2 )m'dm' 2 . (20) 



Here p is the mass density, stratified on surfaces of constant m 2 = w 2 + z 2 /(l — e 2 ), and 
e = J 1 — b 2 /a 2 . Using our adopted b/a = 0.879 for Omega Centauri, and writing S = M 2 S' 



with M the total mass, this becomes W = —31.8GM 2 S'. We can evaluate S' using the 
number density profile derived in §3; the result is S' = 1.463 x 10~ 3 pc^ 1 , assuming a 
distance to Omega Centauri of 5.2 kpc. 

The total kinetic energy is K = 2K X + K z , where (as above) the z- axis is perpendicular 
to the equatorial plane. The tensor virial theorem gives 

K * = (%r) K * ( 21 ) 



\w xx 



with 

W zz (b\ 2 



0.9010 (22) 



( IRoberts 1961 ). Thus K = 2.901^ = 1.451M(V^), with (V?) the (mass weighted) mean 
square line-of-sight velocity of the stars. Equating 2K and \W\ then gives 

M = 1.445 x 1O 4 (V^)M , (23) 

where (V 2 ) is measured in km s _1 . This estimate of the mass is independent of any 
assumptions about the internal kinematics, though it is strongly dependent on the 
assumption that mass follows light. 
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Fig. L— 

Positions of stars with velocities measured by CORAVEL in uj Centauri. Distances are 
measured in minutes of arc from the cluster center. Filled/open circles are stars with radial 
velocities that are greater/less than the cluster mean. The solid line shows the orientation 
of the photometric minor axis as determined by White & Shawl (1987); this is the Z axis 
in what follows. The plotted ellipse has a mean radius equal to the uj Centauri core radius 
as determined by Peterson & King (1975), and an eccentricity equal to the average value 
for uj Centauri, e = 0.121, as determined by Geyer et al. (1983). 

Fig. 2.- 

Surface brightness (a) and space density (b) profiles for uj Centauri. The data in (a) 
are from Table 1 of Meylan (1987); the solid line is the solution to the optimization problem 
defined by Eq. (Q). The dashed lines in (b) are 95% confidence bands on the estimate of 
v{r). Both profiles are normalized to unit total number. 

Fig. 3.- 

(a) The line-of-sight, rotational velocity field of uj Centauri, obtained as the solution to 
the optimization problem (|4]). The Z-axis is parallel to the isophotal minor axis and east 
is toward the left; distances are measured in arc minutes. Dots indicate positions of stars 
in the kinematical sample. Contours are labelled in km s _1 , measured with respect to the 
cluster mean; dashed contours indicate negative velocities, (b) The line-of-sight velocity 
dispersion field of uj Centauri, obtained as the solution to the optimization problem QT2D. 
Plotted are contours of the velocity dispersion about the local mean shown in (a); velocity 
errors have also been removed in quadrature, (c) Heavy line: line-of-sight rotational velocity 
along the major axis of uj Centauri, in km s _1 ; thin line: antisymmetrized profile, (d) 
Heavy line: line-of-sight velocity dispersion along the major axis of uj Centauri; thin line: 
symmetrized profile. 

Fig. 4.— 

Like Fig. 3, except that the perspective rotation due to the proper motion of uj 
Centauri has been removed from the measured velocities. 

Fig. 5.- 

Estimate of the mean azimuthal velocity in the meridional plane of uj Centauri, 
obtained as the solution to the optimization problem (§). Distances are in arc minutes and 
contours are labelled in km s _1 . 
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Fig. 6.— 

(a) Mean azimuthal velocity in the equatorial plane, b) Mean azimuthal velocity in 
the meridional plane as a function of z at w = 7'. Dotted lines are 95% confidence bands, 
computed via the bootstrap. Distances are expressed in minutes of arc. 

Fig. 7.- 

Estimates of the meridional (a) and azimuthal (b) velocity dispersions, a and a^, in 
the meridional plane, obtained as solutions to the constrained optimization problem defined 
by equations ([T2|) and fll3|). Distances are expressed in minutes of arc and contours are 
labelled in km s _1 . 

Fig. 8.— 

Estimates of the velocity dispersion profiles along the major (a) and minor (b) axes in 
the meridional plane. Thick lines: a\ thin lines: o$. Dotted lines are 95% confidence bands 
computed via the bootstrap. Distances are expressed in minutes of arc. 

Fig. 9.— 

Estimates of the velocity anisotropy 1 — cr^/cr 2 along the major (a) and minor (b) 
axes in the meridional plane. Dotted lines are 95% confidence bands computed via the 
bootstrap. Distances are expressed in minutes of arc. 

Fig. 10 — 

Estimates of the gravitational potential (a) and the mass density (b) derived from the 
Jeans and Poisson equations. Thin lines are the corresponding functions derived under the 
assumption that mass follows light in u Centauri. Distances are expressed in minutes of arc. 

Fig. 11 — 

Dependence of $ and p on distance along the major axis in uo Centauri. Solid lines 
are the dynamically-derived estimates, with 95% confidence bands (dashed). Thin lines 
are the corresponding functions derived under the assumption that mass follows light in uj 
Centauri. Units of p are M Q pc~ 3 ; the potential is expressed in units such that the total 
mass of u Centauri is unity. Distances are expressed in minutes of arc. 

Fig. 12 — 

Two-integral distribution functions for uj Centauri. (a) f + (E,Lz); (b) f-(E,L z ). 
Heavy lines are the curves of maximum L z as a function of E. Contours are separated by 
0.27 in log 10 /. E and L z are expressed in units such that the total mass of u Centauri is 
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Fig. 2.- 
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